!
! tester - astrometry
!
!
! $ gfortran -Wall -g -p -fbounds-check -fimplicit-none projections.f95 atrafo.f95 \
!            lsfit.f95 tester.f95 -L../lib -lminpack -lm
!
! Copyright © 2011 F.Hroch (hroch@physics.muni.cz)
!
! This file is part of Munipack.
!
! Munipack is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
! 
! Munipack is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with Munipack.  If not, see <http://www.gnu.org/licenses/>.
!

program tester

  use astrotrafo
  use lsfit
  use absfit
  use robfit

  implicit none

  character(len=*), parameter :: type = 'GNOMONIC'
  integer, parameter :: ns = 9
  integer, parameter :: dbl = selected_real_kind(15)
  real(dbl), dimension(ns) :: x,y,a,d,a2,d2,x0,y0
  real(dbl), parameter :: a0 = 330.0, d0 = 30.0, r = 0.6
  real(dbl), parameter :: scale0 = 1.8e3, pa0 = 0.0
  real(dbl), parameter :: w = 1000, h = 800
  real, parameter :: ge = (2.0/3600.0)*scale0

  real(dbl) :: q,acen,dacen,dcen,ddcen,sc,dsc,pa,dpa,rms,s0,sres,mad
  integer :: i
  type(TypeProjection) :: t

  do i = 1, ns
     call random_number(q)
     a(i) = a0 + r*(q - 0.5) 
     call random_number(q)
     d(i) = d0 + r*(q - 0.5)
  end do

  t%type = type
  t%acen = a0
  t%dcen = d0
  t%xcen = w/2.0
  t%ycen = h/2.0
  t%scale = scale0
  t%pa = pa0

  call trafo(t,a,d,x0,y0)

  do i = 1, ns
     x(i) = x0(i) + gnoise(ge)
     y(i) = y0(i) + gnoise(ge)
!     write(*,*) gnoise(0.14),gnoise(0.14)
  end do

  do i = 1,ns
     write(*,'(2f10.5,2f20.3,2f10.5)') a(i),d(i),x(i),y(i)
  end do

  acen = a0+0.1
  dcen = d0-0.1
  sc = 1.0
  pa = 0
!  call lsmin(type,a,d,w/2.0,h/2.0,w,h,x,y,acen,dacen,dcen,ddcen,sc,dsc,pa,dpa,rms,s0)
  call absmin(type,a,d,w/2.0,h/2.0,w,h,x,y,acen,dcen,sc,pa,rms,s0)
  mad = rms
  call robmin(type,a,d,w/2.0,h/2.0,w,h,x,y,mad,acen,dacen,dcen,ddcen,sc,dsc,pa,dpa,rms,s0)

  write(*,99) acen,dacen,dacen/acen
  write(*,99) dcen,ddcen,ddcen/dcen
  write(*,99) sc,dsc,dsc/sc
  write(*,99) pa,dpa,dpa/pa
  write(*,*) rms,rms*3600
  write(*,*) s0,s0*3600
99 format(f20.10,'+-',f20.10,'  ',1es20.5)

  ! residuals
  t%acen = acen
  t%dcen = dcen
  t%scale = sc
  t%pa = pa
  sres = 0.0
  call invtrafo(t,x,y,a2,d2)
  do i = 1,ns
     write(*,'(2f10.5,2f10.5,2f8.3)') a(i),d(i),a2(i),d2(i), &
          3600*(a(i)-a2(i))/cos(dcen/57.3),3600*(d(i)-d2(i))
     sres = sres + ((a(i) - a2(i))/cos(dcen/57.3))**2 + (d(i) - d2(i))**2
  end do

  write(*,*) sqrt(sres/ns/2.0)*3600.0

contains

   function gnoise(noise)

    real :: gnoise, noise
    real :: x

    call random_number(x)
    gnoise = noise*invdist(x)
    

  end function gnoise


  function invdist(xx)

    real :: invdist,xx

    ! rational approximation of an inverse to a cumulative function
    ! of Gaussian distribution  with precision better than 0.00045
    ! J.Andel: Statistical methods, Matfyz Press, Prague 1991
    
    real :: w,f,x
    logical :: interval

    x = xx
    if( x < 0.0 ) then
       invdist = 0.0
    elseif( x > 1.0 )then
       invdist = 1.0
    else
       interval = x < 0.5
       if( .not. interval ) x = 1.0 - x + epsilon(1.0)
       w = sqrt(-2.0*log(x));
       f = -w + (2.515517 + w*(0.802853 + w*0.010328))/ &
            (1.0 + w*(1.432788 + w*(0.189269 + w*0.001308)));
       if( interval ) then
          invdist = f
       else
          invdist = -f
       endif
    endif
    
  end function invdist


end program tester
